#!/usr/bin/perl
###############################################################################
#
# Provides an algorithm for resistance antiviral calculus
# Developed by Eduardo Sanchez [eduardo.sanchez@fudepan.org.ar]
#
###############################################################################
#
# FuDePAN
# Copyright (C) Eduardo Sanchez, 2008-2009, FuDePan
#
#   min-retro-calc is free software: you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
#   min-retro-calc is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with min-retro-calc.  If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################

use strict;
use warnings;
use Getopt::Long;
require 5;

# some definitions
my $splitted_token="_SPLITTED_";

use constant ZFALSE => 0;
use constant ZTRUE => 1;

# unbuffer the outputs
select STDERR; $| = 1;
select STDOUT; $| = 1;

#Verify the following table.. it was taken from the original rnaCodonTable
my %rnaCodonInverseTable=
(
	"A" => ["gcg","gct","gca","gcc"],
	"R" => ["agg","cgc","cgt","cga","aga","cgg"],
	"N" => ["aac","aat"],
	"D" => ["gat","gac"],
	"C" => ["tgt","tgc"],
	"Q" => ["caa","cag"],
	"E" => ["gaa","gag"],
	"G" => ["ggc","ggg","ggt","gga"],
	"H" => ["cat","cac"],
	"I" => ["att","ata","atc"],
	"L" => ["ctt","ctg","tta","cta","ctc","ttg"],
	"K" => ["aaa","aag"],
	"M" => ["atg","aug"],
	"F" => ["ttc","ttt"],
	"P" => ["cca","ccc","ccg","cct"],
	"S" => ["tct","tcc","tca","tcg","agt","agc"],
	"T" => ["acg","aca","acc","act"],
	"W" => ["tgg"],
	"Y" => ["tac","tat"],
	"V" => ["gta","gtc","gtg","gtt"]
);

my %acceptedNucleotids=
(
  "a" =>[],
  "c" =>[],
  "g" =>[],
  "t" =>[],
  "u" =>[]
);

my $threshold = 60;
my $verbose;
my $callerScriptName = $0; #Grabbing out the caller script name

&main();

# Main Algorithm Section

# -----------------------------------------------------------
#            Triplet Handling Section
# -----------------------------------------------------------

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: dist_triplet
# Arguments: 1) String (length 3) of 1st triplet
#            2) String (length 3) of 2nd triplet
# Return: the distance between the 1st and 2nd triplet (if equal returns 0)
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub dist_triplet($$){
 my $distance=3;
 my $triplet_1 = shift(@_);
 my $triplet_2 = shift(@_);
 if ((length($triplet_1)== length($triplet_2)) and (length($triplet_1)==3))
  {
    for (my $index=3;$index>0;$index--)
       {
          if (chop($triplet_1) eq chop($triplet_2)) { $distance--; }
       }
  } else
  {
     die "Error: please verify the length arguments of dist_triplet function.";
  }

 return $distance;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: getTripletsOfAA
# Argument: String coding the amino acid (for instance: "A","Y", etc).
# Return: Array with all the triplets of a particular amino acid.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub getTripletsOfAA($){
    my $amino_acid= shift(@_);
    my @result=();
    #Verify if the amino acid has any triplet associated.
    if (defined(@{ $rnaCodonInverseTable{$amino_acid} }))
      {
        @result = @{ $rnaCodonInverseTable{$amino_acid} };
      }
    return @result;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: mindist
# Arguments: 1) String coding a triplet
#            2) String coding an Amino Acid (i.e "A")
#            3) Ocurrencies Accumulator
# Return: lower distance between the triplet and all the triplets coding
#         the selected amino acid.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub mindist($$$){
  my $triplet=shift(@_);
  my $amino_acid=shift(@_);
  my $accumulator=shift(@_);
  my @set_of_triplet = getTripletsOfAA($amino_acid);
  my $distance=3;
  my $aux_distance;
  my $acc=0;

  fudepan_print("triplets of aa $amino_acid -> @set_of_triplet\n");

  for (my $index=0;$index<=$#set_of_triplet;$index++)
     {
        $aux_distance = dist_triplet($triplet,$set_of_triplet[$index]);
        if ($distance == $aux_distance)
        { #print "aux_distance = distance\n";
          $acc++;
        }  elsif ($distance > $aux_distance)
        {
          #print "set aux_distance=$aux_distance instead of $distance\n";
          $distance = $aux_distance;
          $acc=1; # one element with that distance..
        }
     }

  $$accumulator = $acc;
  return $distance;
}


# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: check4ValidNucleotids
# Arguments: 1) string with sequence of nucleotids (or just a triplet will be fine)
# Return: ZTRUE if it's ok, or ZFALSE otherwise
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *


sub check4ValidNucleotids($)
{
  my $string = lc shift(@_);
  my $char;
  my $result = ZTRUE;
  while (($string ne "") and $result)
  {
    $char = chop($string);
    if (not (exists($acceptedNucleotids{$char})))
     {
       $result = ZFALSE;
     }
  }
  return $result;
}


# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: triplete
# Arguments: 1) Result string reference
#            2) Sequence string
#            3) First element index
#            4) Hash of dash positions
#            5) Nucleotid relationship tables (original and updates)
#            6) Antiviral
# Return: ZTRUE if it's ok, or ZFALSE otherwise
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub triplete
{
   my ($triplet,$sequence,$position,$dashPositionsRef,$NT_NTUpdatePosRef,$antiviral)= @_;
   my $result = ZFALSE;
   my @offset_chars;
   my $numberOfoffset_chars=0;
   my $tripletSafePosition;
   my @safePositions=();
   my $base_position = $position;
   my %hash_keys = map { $_ => $_ } @{$dashPositionsRef}; #convert the array to hash in order to use exists function

   @offset_chars = grep { $_ < $position } @{$dashPositionsRef};
   # the position start will be.. some position + $# offset_chars..
   # offset_chars are the used character by the previuos triplets
   if (@offset_chars)
   {
     $numberOfoffset_chars = $#offset_chars+1;
   }

   while($numberOfoffset_chars>0)
   {
     while (exists($hash_keys{$position}))
     {  $position++;
     }
     $numberOfoffset_chars--;
     $position++;
   }
   # until here, I "eat" safe positions already used by previous triplets..
   # i have to check if following positions are safe or not..

   for ($tripletSafePosition=0;$tripletSafePosition<3;$tripletSafePosition++)
     {  #update zero safe position..
       while (exists($hash_keys{$position}))
       {  $position++;
       }
      # print "position $tripletSafePosition safe: $position\n";
       $safePositions[$tripletSafePosition]= $position;
       $position++;
     }

    $NT_NTUpdatePosRef->{$antiviral}->{$base_position}=$safePositions[0];

   #now, check if safePos are inside the boundaries.. (just the lastone..)
   if ($safePositions[2] < length ($sequence))
   {
   #since positions are inside the sequence.. get the triplet..
    $$triplet = "";
    for ($tripletSafePosition=0;$tripletSafePosition<3;$tripletSafePosition++)
    {
      $$triplet = $$triplet . substr($sequence,$safePositions[$tripletSafePosition],1);
    }
    $result = check4ValidNucleotids ($$triplet);
    #fudepan_print("triplete ->$$triplet<-\n");
   }

   return $result;
}


# -----------------------------------------------------------
#                         Main Section
# -----------------------------------------------------------

sub main()
{
  my $NTposition;
  my $antiviral;
  my %resultAccMatrix = ();
  my %sequenceTable = ();
  my %ResistanceTable = ();

  my $AA;
  my $sequence;
 # Auxiliary Variables..
  my $singleQuoteP; # P'
  my $auxiliaryTriplet= "";
  my $auxiliaryMinDist;
  my $auxiliaryMin;
  my $aux_check;
  my @sequenceIdTable; # contains the sequences following the sequences file order
  my $sequenceId;
  my $helpMode;
  my $merge;
  my $noSplit;
  my $autoExpand;
  my $firstIsMaster;
  my $excelFile;
  my @resistanceTableFile;
  my $sequencesFastaFile;
  my $masterSequence="";
  my $accumulatedDistance;
  my %NT_NTUpdatePos;
  my $displayPositions;
  my @dashPositions=();
  my %results=();

 #Parsing out the provided arguments
 #Check that enough arguments were provided

 if ( @ARGV > 0 )
{
	#Display help if one of the needed arguments are missing
	GetOptions(
            'first-is-master' => \$firstIsMaster,
            'excel-file:s' => \$excelFile,
            'resistance-table-file=s@' => \@resistanceTableFile,
            'sequences-fasta-file=s' => \$sequencesFastaFile,
            'display-positions:s' => \$displayPositions,
            'merge' => \$merge,
            'no-split' => \$noSplit,
            'auto-expand' => \$autoExpand,
            'verbose' => \$verbose,
	    'help|?' => \$helpMode,

	    ) or usage($callerScriptName, "Invalid Arguments were provided");
}
else
{
	usage($callerScriptName, "No arguments were provided");
}

# --first-is-master --outputFile edu.txt --resistance-table-file Scores_EDU.txt --sequences-fasta-file fasta.txt

#Help menu was called
if( defined ($helpMode))
{
	usage($callerScriptName);
}

# Required Arguments
  if (not (@resistanceTableFile))
  {
   die "You must provide at least one resistance table file using option \"--resistance-table-file\".\n";
  }

  if (not (defined $sequencesFastaFile))
  {
    die "You must provide one sequence fasta file using option \"--sequences-fasta-file\".\n";
  }


#Optional Arguments

if (defined ($displayPositions))
{
  if (not(($displayPositions eq "AA") or ($displayPositions eq "NT")))
  {
   $displayPositions="BT"
  }
}
else
{
  $displayPositions="BT"; #BT = Both, aminoacid and nucleotid position.
}



  if (not (defined $firstIsMaster))
  {
     $firstIsMaster=0;
  }

  # Initialize Resistance Table
  # NEVER reorder populateResistanceTable, createRangesFile or get SequencesFromFastaFile!
  populateResistanceTable(\%ResistanceTable,@resistanceTableFile);

  if ( defined($autoExpand))
   {
     createRangesFile(\%ResistanceTable);
     my $result = `fasta-expander.py -f "$sequencesFastaFile" -g -r auto_ranges.tmp > expanded.tmp`;
     $sequencesFastaFile = "expanded.tmp"; 
   }

  # Initialize Hash Table of Sequences
  getSequencesFromFastaFile(\@sequenceIdTable, \%sequenceTable, $sequencesFastaFile, $firstIsMaster, \$masterSequence);

  preprocess_resistance_table_with_master_sequence(\%ResistanceTable,$masterSequence,\@dashPositions,$firstIsMaster);


  fudepan_print("Aminoacid & Nucleotid Resistance Table\n");
  showResistanceTable(\%ResistanceTable,"BT",\%NT_NTUpdatePos);


  # end initialization
  print "Running $callerScriptName script...\n";
  %results=();
 for (my $index=0;$index<=$#sequenceIdTable;$index++)
 {

  $sequenceId = $sequenceIdTable[$index];
  $sequence = $sequenceTable{$sequenceId};

   foreach $antiviral ( keys %ResistanceTable ) #foreach antiviral 3TC, AZT..
   {
     foreach $NTposition ( sort by_value keys %{$ResistanceTable{$antiviral}} ) #for each antiviral's position P
     {
        $results{$antiviral}{$sequenceId}{$NTposition}[0]=3; # algorithm results
        $results{$antiviral}{$sequenceId}{$NTposition}[1]=0; # occurrencies number

        $resultAccMatrix{$NTposition}{$antiviral}[0]=3; # algorithm results
        $resultAccMatrix{$NTposition}{$antiviral}[1]=0; # occurrencies number

        $aux_check = triplete(\$auxiliaryTriplet,$sequence,$NTposition,\@dashPositions,\%NT_NTUpdatePos,$antiviral);
        if ($aux_check == ZFALSE)
          { my $fixedPosition = $NTposition+1;
            print "Error in sequence $sequenceId: triplet \"$auxiliaryTriplet\" at position: $fixedPosition.\n";
            $results{$antiviral}{$sequenceId}{$NTposition}[0]=-1;
          }
          else
          {
            foreach $AA ( @{$ResistanceTable{$antiviral}{$NTposition}})
            {
              $auxiliaryMinDist=mindist($auxiliaryTriplet,$AA,\$accumulatedDistance);
              if ($auxiliaryMinDist < $results{$antiviral}{$sequenceId}{$NTposition}[0])
              {
               $results{$antiviral}{$sequenceId}{$NTposition}[0]=$auxiliaryMinDist;
               $results{$antiviral}{$sequenceId}{$NTposition}[1]=$accumulatedDistance;
              } elsif ($auxiliaryMinDist == $results{$antiviral}{$sequenceId}{$NTposition}[0])
                 {
                   $results{$antiviral}{$sequenceId}{$NTposition}[1] = $results{$antiviral}{$sequenceId}{$NTposition}[1] + $accumulatedDistance;
                 }
            }
          }
     }

   }


   showResistanceTable(\%ResistanceTable,$displayPositions,\%NT_NTUpdatePos);

 }


  if ( defined($autoExpand))
     {
       my $removeSplit="yes";
       mergeSplittedSequences(\@sequenceIdTable,\%results,$removeSplit); #merge all results from SPLITTED sequecies
       my $result = unlink("expanded.tmp"); #remove temporary files
       $result = unlink("auto_ranges.tmp");
     }
  else  # Begin changes for merging results...
  {
    if( defined ($merge))
    {
      mergeSplittedSequences(\@sequenceIdTable,\%results,$noSplit);
    }

  } 

  showResults(\%results,\@sequenceIdTable,\%NT_NTUpdatePos,$displayPositions,$excelFile);
 # end changes for merging results...



  print "$callerScriptName script finished.\n";

}

# -----------------------------------------------------------
#   Merge section
# -----------------------------------------------------------

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: mergeSplittedSequences
# Arguments: 1) Sequence id table
#            2) Master Results table
# Return: both entries updated.-
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub mergeSplittedSequences(\@\%$)
{
   my $sequenceIdTableRef = shift (@_);
   my $resultsRef = shift (@_);
   my $noSplit = shift (@_);
   my %splittedSequences = ();
   my $splitted_index;
   my $sequence_name;

   # check if "SPLIT" tag is inside the sequence_name
   for (my $index=0;$index< $#{ $sequenceIdTableRef }+1;$index++)
   {
     $sequence_name = $sequenceIdTableRef->[$index];
     $splitted_index= index($sequence_name,$splitted_token);
     if ($splitted_index > 0)
     {  my $total_seq_name_length = length $sequence_name;
        my $original_sequence_name = substr($sequence_name,0,$splitted_index);
        my $additional_data = substr($sequence_name,$splitted_index+length $splitted_token ,$total_seq_name_length);
        $splittedSequences{$original_sequence_name}{$sequence_name} = $additional_data;
     }
   }

  # at this point, in the splittedSequences hash, we have the original_sequence and the rest of his splitted names
  foreach my $seq_original_name (sort keys %splittedSequences)
  {
     foreach my $seq_split_name (keys %{$splittedSequences{$seq_original_name}})
     {
       #here I know the splitted sequence name..
       # we have to retrieve the old splitted values and store the newones if neccessary
       foreach my $antiviral (keys %{$resultsRef})
       {
         foreach my $NTposition (keys %{$resultsRef->{$antiviral}->{$seq_split_name}})
         {
           my $temp_result = $resultsRef->{$antiviral}->{$seq_split_name}->{$NTposition}->[0];
           my $temp_occurrencies = $resultsRef->{$antiviral}->{$seq_split_name}->{$NTposition}->[1];

           if (not(exists $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}))
           {
             $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[0] = $temp_result;
             $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[1] = $temp_occurrencies;
           } else
           { # there is an entry for this splitted sequence..
            my $actual_result = $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[0];
            my $actual_occurrencies = $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[1];
            if ($temp_result < $actual_result)
            {
              $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[0] = $temp_result;
              $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[1] = $temp_occurrencies;
            } elsif (($temp_result == $actual_result) and ($temp_occurrencies > $actual_occurrencies))
            {
              $resultsRef->{$antiviral}->{$seq_original_name}->{$NTposition}->[1] = $temp_occurrencies;
            }
           } #close else
         } # close foreach NT
       } #close foreach antiviral
     } #close foreach $seq_split_name

    # push @{$sequenceIdTableRef}, $seq_original_name; #add the original sequence to results table
     # if needed, here is the place to remove the _SPLITTED_ sequences in sequenceIdTableRef.
 } #close foreach $seq_original_name

 reorderSplittedSequences($sequenceIdTableRef,$noSplit);

}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: reorderSplittedSequences
# Arguments: 1) Sequence id table
#            2) noSplit flag
# Return: the sequence id table will be updated with the original name
#         just before the first splitted name.
#         in case that noSplit flag active, remove the splitted sequence name
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub reorderSplittedSequences(\@$)
{  my $sequenceIdTableRef = shift (@_);
   my $noSplit = shift (@_);
   my %originalSplittedNames = ();
   my @sequenceIdTableUpdated = ();
   my $original_sequence_name;
   # check if "SPLIT" tag is inside the sequence_name
   for (my $index=0;$index< $#{ $sequenceIdTableRef }+1;$index++)
   {
     my $sequence_name = $sequenceIdTableRef->[$index];
     my $splitted_index= index($sequence_name,$splitted_token);
     if ($splitted_index<0)
     { # there is no issues in adding the sequence name
       push @sequenceIdTableUpdated, $sequence_name;
     } else
     {
       $original_sequence_name =  substr($sequence_name,0,$splitted_index);
       if (not(exists $originalSplittedNames{$original_sequence_name}))
       { # add the sequence name to the existing list..
         $originalSplittedNames{$original_sequence_name}=$original_sequence_name;
         # push the sequence name.. just once..
         push @sequenceIdTableUpdated, $original_sequence_name;
       }
       if (not(defined $noSplit))
       { #only push the sequence splitted name, if the flag is not true
         push @sequenceIdTableUpdated, $sequence_name;
       }
     }
   }
 #replace the global sequence id table with the local table
 @{$sequenceIdTableRef} = @sequenceIdTableUpdated;

}

# -----------------------------------------------------------
#                         Auxiliary Section
# -----------------------------------------------------------


# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: preprocess_resistance_table
# Arguments: 1) Resistance Table Hash Reference
#            2) Master sequence string
#            3) hash of dashes position
#            4) first-is-master flag
# Return: * hash of dashes table populated.
#         * resistance table position in nucleotids.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub preprocess_resistance_table_with_master_sequence(\%$\@$)
{ my $resistanceTableRef = shift @_;
  my $masterSequence = shift @_;
  my $dashPositionsRef = shift @_;
  my $firstIsMaster = shift @_;
  my %resistanceTableUpdated = ();
  my $NTposition;
  my $AAposition;
  my $where=0;

  if ($firstIsMaster)
  {  fudepan_print("Resistance Table will be updated with offsets...(first-is-master flag)\n");
     # find all the "-" in sequence
     while ($where>-1)
     {  $where = index($masterSequence, "-",$where);
        if ($where>-1)
        {   push(@{$dashPositionsRef},$where);
            $where++;
        }
     }
  }

 # convert just Aminoacid position to NT position.. (whitout taking care "-")
  foreach my $antiviral ( keys %{$resistanceTableRef} )
  {
   # print "antiviral $antiviral\n";
    foreach my $position ( sort by_value keys %{$resistanceTableRef->{$antiviral}} ) #para cada posicion P del Antiviral
    {
     $AAposition= $position;
     $NTposition = aminoacid2NucleotidPosition($position); # nucleotid position.. (0 to ..)
    # print "NT conversion: $NTposition.. \n";
     $resistanceTableUpdated{$antiviral}{$NTposition} = $resistanceTableRef->{$antiviral}->{$AAposition};
    }
  }

   # update the original resistance table..
   %{$resistanceTableRef} = %resistanceTableUpdated;

 #here, the ResistanceTable was successfully updated with the master relative position.
 return 0;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: getSequencesFromFastaFile
# Argument: 1) Sequence Id hash reference (for output)
#           2) Sequence Table hash Reference (hash table with  "sequence Id" -> "sequence")
#           3) fasta file name (input)
#           4) is_first_master flag (for input)
#           5) masterSequence (for output)
# Precondition: well-formed fasta File:
#               >sequence name <enter>
#               1 or more lines with sequence content
# Optional: first sequence is master:  >sequence name <enter>
#                                      1 master sequence
#                                      rest of sequences..
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub getSequencesFromFastaFile
{   my $sequenceIdRef = shift @_; #since the hash has no order, we need to take care of him in the future
    my $sequenceTableRef = shift @_;
    my $fastaFileName = shift (@_);
    my $is_first_master = shift (@_);
    my $masterSequence = shift (@_);
    my $seqId="noSeqId";

    open FASTA_FILE, "< $fastaFileName" or die "Could not open $fastaFileName\n";
    my @fileContents = <FASTA_FILE>; # put in memory the complete file content
    close FASTA_FILE;
    chomp(@fileContents);
    @{$sequenceIdRef} = ();

 # get the master sequence if the flag is activated...
    if ($is_first_master)
    {
       shift @fileContents; # remove label
       $$masterSequence = lc shift @fileContents;
    }

    #    parse the rest of the fasta file...
    foreach my $aux_line (@fileContents)
     {
       if ($aux_line =~ m/>/)
        {  # remove ">" from sequence id
           $seqId= substr($aux_line,1);
           push (@{$sequenceIdRef}, $seqId);
        }
      else
       { if (defined($sequenceTableRef->{$seqId}))
          {
            $sequenceTableRef->{$seqId} = $sequenceTableRef->{$seqId} . lc ($aux_line);
             # For multi-lines sequences, convert to lowercase
          } else
          {
            $sequenceTableRef->{$seqId} = lc ($aux_line);
          }
       }
    }

}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: showResultMatrix
# Argument: 1) hash of matrix result
#           2) nucleotid relationship table
#           3) display type (NT, AA, BT)
#           4) sequenceId
#Results: Display through standard output the algorithm results.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub showResults(\%\@\%$$)
{
   my ($matrix_ref, $sequenceIdTable, $NT_NTUpdatePosRef,$displayType,$excelFile) = @_;
   my $sequenceId;
   my $nt_updated_pos;
   my $aaPosition;
   my $display_header=0;
   my $positionsLine;
   my $resultsLine;

   my $label;
   if ($displayType eq "NT")
   {
    $label = "nucleotid position";
   } elsif ($displayType eq "AA")
   {
    $label = "aminoacid position";
   } else
   {
    $label = "[aa pos,nt pos]";
   }

   print ("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n\n");

   if (not(%$matrix_ref))
   {
     print "The result table is empty.\n";
   }
   else
   {
     if (defined ($excelFile))
     {
       open TAB, ">$excelFile.csv" or die "Couldn't open file. $!\n";
     }
     foreach my $antivir ( sort { $matrix_ref->{$a} <=> $matrix_ref->{$b}} keys %$matrix_ref )
     {
       print "$antivir\n";
      if (defined ($excelFile))
      {
        print TAB "Antiviral: $antivir\n";
      }

       $display_header=1;
       for (my $index=0;$index<=$#{$sequenceIdTable};$index++)
       {
         $sequenceId = $sequenceIdTable->[$index];
         $positionsLine="";
         $resultsLine="";
         foreach my $pos ( sort by_value keys %{$matrix_ref->{$antivir}->{$sequenceId}} )
         {
          $nt_updated_pos = $pos;
          $aaPosition = nucleotid2AminoacidPosition($pos);
          if (exists ($NT_NTUpdatePosRef->{$antivir}->{$pos}))
            { $nt_updated_pos=$NT_NTUpdatePosRef->{$antivir}->{$pos}; #if the nt position was updated, display the last one
            }
            $nt_updated_pos++; # add one in order to start position at 1.

            if ($displayType eq "NT")
            {
              $positionsLine = $positionsLine . "\t" . "$nt_updated_pos";
            } elsif ($displayType eq "AA")
            {
              $positionsLine = $positionsLine . "\t" . "$aaPosition";
            } else
            {
              $positionsLine = $positionsLine . "\t" . "[$aaPosition, $nt_updated_pos]";
            }


            if (defined ($matrix_ref->{$antivir}->{$sequenceId}->{$pos}))
             {
               $resultsLine=  $resultsLine . "\t" . "$matrix_ref->{$antivir}->{$sequenceId}->{$pos}->[0]($matrix_ref->{$antivir}->{$sequenceId}->{$pos}->[1])";
             }
         }

         $positionsLine = "sequence/$label: $positionsLine";
         $resultsLine = $sequenceId . $resultsLine;
         if ($display_header)
         {
           print "$positionsLine\n";
           if (defined ($excelFile))
           {
             print TAB "$positionsLine\n";
           }
           $display_header=0;
         }
         print ">$resultsLine\n";
         if (defined ($excelFile))
         {
           print TAB ">$resultsLine\n";
         }
       }
     }
   if (defined ($excelFile))
   {
     close TAB;
   }
  } #else

}


#begin sorting hash functions section

sub by_value {$a <=> $b}

#end sorting hash functions section

# -----------------------------------------------------------
# Begin Read Table Section
# -----------------------------------------------------------

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: populateResistanceTable
# Arguments: 1) resistance table hash
#            2) List of Table File Names...
# Purpose: populate the Resistance Table for future handle.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub populateResistanceTable
{
  my ($resistanceTableRef, @file_table_names) = @_;
  if ($#file_table_names<0)
  {
    die "You must provide at least one file with the antiviral resistance.\n";
  }
  foreach my $file_name (@file_table_names)
  {
     createOrModify_ResistanceTableFromFile($resistanceTableRef,$file_name);
  }
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: createOrModify_ResistanceTableFromFile
# Argument: 1) resistance table reference
#           2) Table File Name...
# Purpose:
# Read file with Antiviral Resistance Table Source, inserting data in
# RAM Resistance Table (only values upper bound threshold).
# This function is the core function of the file.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub createOrModify_ResistanceTableFromFile($$){
    my $resistanceTableRef = shift (@_);
    my $table_filename = shift(@_);

   # print "Openning file $table_filename...\n";
    open TABLE_FILE, "< $table_filename" or die "Could not open $table_filename\n";

    # Read Header Table ...
    my $header_line = <TABLE_FILE>;
    chomp($header_line); #remove chr(13);

    my @headers_arr = split /\s+/, $header_line; # all headers space-separated.
    my $master_table_length = $#headers_arr; # master table length

    # Read Body Content Table...

    my @fileContents = <TABLE_FILE>;
    chomp(@fileContents);
    foreach my $aux_line (@fileContents)
    {
       if (defined($aux_line) and ($aux_line ne ""))
       {
         my @aux_arr = split /\s+/, $aux_line;
         if ($master_table_length != $#aux_arr)
         {
           die "File $table_filename has not equal column length. Header = $master_table_length <> Current = $#aux_arr.\n";
         }

         for (my $aux_index=2; $aux_index<$#aux_arr+1; $aux_index++)
         {
          if ($aux_arr[$aux_index] >= $threshold)
          {
            # Position = $aux_arr[0];
            # AminoAcid = $aux_arr[1];
            # Antiviral = $headers_arr[$aux_index];
            # Local Threshold Value = $aux_arr[$aux_index]
            my $aminoacid_name = $aux_arr[1];

            if (checkForWrongAminoAcidName($aminoacid_name) == ZTRUE) # check for "-" or "insert"
            {
               insertResistanceTable($resistanceTableRef,$headers_arr[$aux_index],$aux_arr[0],$aminoacid_name);
            }
            else
            {
               fudepan_print("WARNING: The file $table_filename contains an invalid amino acid name: $aminoacid_name! \n");
            }
          }
         }
       }
    }
   close TABLE_FILE;

}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
# Auxiliary functions:
# Function name: createRangesFile
# Arguments: 1) Resistance Table hash
# Purpose: Create proper ranges file only for resistance table positions
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *


sub createRangesFile
{
    my $resistanceTableRef = shift (@_);
    my $antiviral;
    my $position;
    my %positions=();
    my $NTposition;


    foreach my $antiviral ( keys %{$resistanceTableRef} )
    {
      foreach my $position ( sort by_value keys %{$resistanceTableRef->{$antiviral}} ) #para cada posicion P del Antiviral
      {
          $NTposition = aminoacid2NucleotidPosition($position); # nucleotid position.. (0 to ..)
          $NTposition++;

          if (not(exists($positions{$NTposition})))
          {
             $positions{$NTposition} = $NTposition; #save unique positions.
          }
       }
    }

    #create ranges file from position Table:
    open RANGES, ">auto_ranges.tmp" or die "Could not create temporary ranges file.\n";
     
    foreach my $disp_pos ( sort by_value keys %positions )
    {
      my $last_disp_pos = $disp_pos+2;
      print RANGES "$disp_pos,$last_disp_pos\n";
    }
    close RANGES;

}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Auxiliary functions:
# Function name: showResistanceTable
# Arguments: 1) Resistance Table hash
#            2) displayType (nucleotid, aminoacid or both)
#            3) position table relationship
# Purpose: Display Antiviral Resistance Table:
#	"3TC" => { 2 => ["A","Y"],
#                   4 => ["K"]
#                 },
#	"AZT" => { 7 => ["A","Y"],
#                   8 => ["K"]
#                 }
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub showResistanceTable
{
    my $resistanceTableRef = shift (@_);
    my $displayType = shift (@_);
    my $NT_NTUpdatePosRef = shift (@_);
    my $antiviral;
    my $position;
    my $index;

    fudepan_print("Resistance Table:\n");

    foreach $antiviral ( sort keys %$resistanceTableRef )
    {
        fudepan_print("$antiviral: { \n");
        foreach $position ( sort by_value keys %{ $resistanceTableRef->{$antiviral} } )
        { my $aaPosition = nucleotid2AminoacidPosition($position);
          my $ntPosition = $position;
          if (exists ($NT_NTUpdatePosRef->{$antiviral}->{$ntPosition}))
          { $ntPosition=$NT_NTUpdatePosRef->{$antiviral}->{$ntPosition}; #if the nt position was updated, display the last one
          }
          $ntPosition++; #display ntPos starting at 1.. not 0
          if ($displayType eq "NT")
          {
             fudepan_print("$ntPosition = ");
          } elsif ($displayType eq "AA")
          {
             fudepan_print("$aaPosition = ");
          } else
          {
            fudepan_print( "($aaPosition, $ntPosition) = ");
          }

          for ($index=0;$index<$#{$resistanceTableRef->{$antiviral}->{$position}};$index++)
          {
            fudepan_print("$resistanceTableRef->{$antiviral}->{$position}[$index],");
          }
           #print the last element without comma.
           fudepan_print("$resistanceTableRef->{$antiviral}->{$position}[$index]\n");
        }
        fudepan_print("}\n");
   }
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: insertResistanceTable
# Arguments: 1) resistance table hash
#            2) Antiviral Name
#            3) Position
#            4) Aminoacid
# Purpose: insert a complete 3-tuple (antiviral, position, aminoacid) in
# memory Resistance Table.
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub insertResistanceTable(\%$$$){
  my $resistanceTableRef = shift(@_);
  my $antiviral = shift(@_);
  my $position = shift(@_);
  my $aminoAcid = shift(@_);
  #my $NTposition = aminoacid2NucleotidPosition($position);;
  #$positionsRef->{$antiviral}->{$NTposition} = $position; # just for tracking the old positions and the newones
  push @{$resistanceTableRef->{$antiviral}->{$position}} , $aminoAcid;
  #fudepan_print("DEBUG: inserting $aminoAcid in position $position of antiviral $antiviral \n");
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub aminoacid2NucleotidPosition($)
{
  my $updated_position = shift;
  $updated_position = ($updated_position-1)*3;
  return $updated_position;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub nucleotid2AminoacidPosition($)
{
  my $updated_position = shift;
  $updated_position = ($updated_position/3)+1;
  return $updated_position;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub checkForWrongAminoAcidName
{  my $aminoAcidName = uc shift(@_);
   my $result = ZFALSE;
   if (exists($rnaCodonInverseTable{$aminoAcidName}))
     {
       $result = ZTRUE;
     }
   return $result;
}

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Function Name: usage
# Argument: 1) script name
#           2) string to be displayed (Optional)
# Return: nothing
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub usage
{
    my ($caller, $stringToDisplay) = @_;

	if ( defined ($stringToDisplay) )
	{
		print STDERR "Problem detected: $stringToDisplay\n";
	}

	#printing help
	die << "EOF";
Resistance calc script
Developed by Eduardo Sanchez [eduardo.sanchez at fudepan.org.ar]
Usage: $caller [options]
Required options:
--sequences-fasta-file	name of the fasta file that lists all the sequences.
--resistance-table-file	name of the file that contains the resistance table. Could be more than one --resistance-table-file option.

Optional options:
--first-is-master	the first sequence will update positions in the resistance table.
--excel-file		name of the output CSV file (without "csv" extension)
--display-positions     show result positions:
                           available options: "AA" -> display original positions
                                              "NT" -> display nucleotid positions of nucleotid
                                              "BT" -> display both, comma separated: i.e. (AA,NT)
--merge                 after processing, merge results of sequencies "nn"_SPLITTED_"mm" in a new "nn" sequence.
--no-split              remove "nn"_SPLITTED_"mm" sequencies from results.
--auto-expand           This option will automatically:
                            1) Create ranges file
                            2) Apply the fasta-expander script
                            3) Merge the results from SPLITTED sequences
                            4) Remove the individual SPLITTED  sequences
                            5) Remove temporary files.
                        Use this option instead of --merge and --no-split options.
--verbose               display all aditional data
--help                  display this help.	
EOF
}


# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Wrapper print to stdout
# Parameter: input to print
# Returns: nothing
# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

sub fudepan_print
{
    my $string = shift;

    if ( defined($verbose) )
     {
     	# flush all output immediately
     	$| = 1;
     	print STDERR "$string";
     }
}
